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This thesis presents a method for the dynamic system identification and simulation 
model development of a small rotary wing UAV. Using aerodynamic parameterization and 
linear state-space modeling techniques, the Bergen Industrial UAV was modeled for 
computer simulation to analyze its inherent stability and control characteristics. The NPS 
designed JANRAD software was utilized to determine the stability and control derivatives 
used in the simulation model. The identification of the UAV dynamic model will aid in the 
development of closed-loop controllers capable of autonomous UAV control. The fidelity 
of the simulation model was verified by comparing the simulation responses with data 
collected from on-board sensors during test flights. 
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L INTRODUCTION 



A. BACKGROUND 

“Everybody ’s simulation model is guilty until proved innocent. ” 

(Thomas H. Lawrence at the 50 th Annual 
Forum of the AHS, Washington, 1994) 

Computer simulation use by the aerospace engineer is beneficial in preliminary 
design, stability and control analysis, and handling qualities determination for air vehicles. 
Automatic control systems enable the stability and handling qualities of an aircraft to be 
augmented to increase performance. With the addition of onboard sensors tied to a 
feedback control system, autonomous flight is possible, and is widely used today on a 
variety of Unmanned Air Vehicles (UAV) within the DoD. Control design theory requires 
a thorough knowledge of the system to be controlled. When dealing with the control of air 
vehicles, an improper design can lead to catastrophic results; therefore, accuracy in the 
characterization of the system to be controlled will facilitate the design process enabling us 
to “model” the aircraft prior to breaking the earthly bounds. The Naval Postgraduate 
School has conducted extensive research in the design of a variety of feedback controllers 
for the Bluebird and FROG fixed-wing UAVs, in addition to the Archytas, VTOL air 
vehicle. To further research in this area, NPS has obtained a Bergen Industrial rotary-wing 
UAV to conduct similar studies in this field by the use of feedback controllers. 



B. APPROACH TO THE PROBLEM 

In the past, control of UAVs was conducted primarily by a ground-based pilot 
using a standard radio controller (RC) where each task had to be monitored by the ground 
pilot for the duration of the flight. With the advent of the global positioning system (GPS), 
the UAV’s position can be tracked electronically vice relying on video imaging projected 
from the UAV to the ground station. Incorporating onboard sensors, feedback controllers 
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or autopilots can now be used to maintain UAV control without full time pilot inputs. The 
Predator UAV can operate autonomously for up to 24 hours using pre-planned routes or 
near-real-time course changes all using onboard sensors to track aircraft position. The 
Naval Postgraduate School has developed UAV controllers which can be controlled 
autonomously or via voice command using a “wearable” voice recognition PC unit. The 
Bergen Industrial Twin UAV (Figure 1.1) was acquired to continue air vehicle control 
research with the ultimate goal of using video imaging to control the autonomous flight of 
a rotary wing air vehicle. Prior to developing a controller for the UAV for autonomous 




Figure L j Bergen Industrial Twin UA V 



flight, it was desired to obtain a mathematical model capable of simulating the air vehicle’s 
response to controller inputs. This paper outlines the effort of the author to obtain an 
accurate model of the UAV by using aerodynamic parameterization. Additionally, it 
provides a method by which a similar modeling and verification process can be conducted 
on similar scale rotary-wing UAVs. It is assumed the reader has an understanding of 
helicopter aerodynamics, aircraft stability and control, and basic system control theory. 
The references recognized throughout this paper provide an excellent source for 
conducting further research. 
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To accomplish the simulation modeling, the following approach was utilized: 

1) Determine UAV physical parameters 

2) Determine stability derivatives 

3) Develop simulation model 

4) Perform flight testing to verify computer model fidelity. 

Using physical measurements, experimental testing and similarity analysis, or 
table look-up, the UAV can be “parameterized” in order to conduct a performance and 
stability analysis. From this physical parameterization of the air vehicle, we wish to 
develop a simulation model capable of conducting frequency and time-response analysis. 
To perform this analysis, it was necessary to characterize the air vehicle system or “plant”. 
To obtain the stability and control derivatives, the NPS developed Joint Army/Navy 
Rotorcraft Analysis and Design (JANRAD) software program was utilized in the 
determination of modeling parameters and frequency response analysis. Finally, the fidelity 
or accuracy of the simulation model can be verified by conducting flight testing with an 
instrumented UAV to obtain aircraft response data for comparison with the simulation 
data. The flight-testing conducted was restricted to a hover due to mechanical and 
structural problems experienced during the research process. It is the goal of the author 
that, through the verification/modification of the simulation model, this work can be used 
for future autonomous controller design. 



3 





































4 



H. MODELING AND SIMULATION 



A. ANALYSIS THEORY 

A mathematical simulation model of a helicopter’s flight dynamics must include the 
important aerodynamic, structural and other integral dynamic effects that combine to 
influence the aircraft's response to pilot inputs. The flight dynamics of the helicopter is 
described by a complex, non-linear, nine degree-of-ffeedom (DOF) system. Describing 
the helicopter’s flight dynamics is a difficult task, but developing an accurate mathematical 
model presents an even greater challenge. A widely used approximation for describing the 
helicopter’s behavior is the linearized, six DOF system. The linearized model is adequate 
for analyzing small perturbations (± 15°) about a trim condition. For an in-depth 
development of the linearized model, the reader is referred to Reference 1. A general 
overview is presented below. 

B. UAV DESCRIPTION 

Future utilization of the selected air vehicle requires the capability to accommodate 
a sizeable payload. Additionally, the initial purchase cost must not be prohibitive within 
the constraints of the research budget. Past UAV projects have benefited from the 
generosity of excess military platforms, such as the Bluebird and the Frog. The Bergen 
industrial UAV was designed as an Industrial platform capable of handling camera 
equipment for aerial photography. Its advertised payload of 20 pounds was within the 
desired scope for handling the onboard telemetry required for the control research 
envisioned by the thesis advisor. 

The UAV designer uses off-the-shelf components for the fuselage and rotor 
systems. In order to increase the UAV’s payload capacity, the standard single-cylinder 
Zenoah engine was modified into a dual-cylinder engine (hence the name Industrial Twin). 
The increased compression (torque) inherent in a dual-cylinder engine was more than the 
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standard (single-cylinder) pull-starter could handle leading to several flight schedule set- 
backs due in part to a broken starter. The RC helicopter is controlled in the same manner 
as a standard single-rotor helicopter; thrust is controlled by main rotor collective pitch, 
longitudinal pitch and roll are controlled by main rotor cyclic (longitudinal and lateral, 
respectfully), and directional control (and anti-torque) are controlled by the tail rotor. The 
UAV is controlled by an outside pilot via remote control using a standard RC transmitter, 
receiver and electric-control servos. Each of the four control inputs receives a separate 
signal to control the appropriate on-board control servo which in turn actuates the 
appropriate control surface (changes blade pitch either cyclically or collectively). The 
transmitted signal from the Futaba® controller is in the form of a pulse width modulation, 
or PWM signal. 

The similarities of the UAV and conventional, full-size helicopter end when we 
discuss the operation of the main rotor system. Like the full scale helicopter, main rotor 
collective is controlled by changing the rotor pitch on all blades simultaneously, or 
collectively, by the uniform movement of the swashplate, but rotor cyclic control is 
slightly different. The Bergen industrial incorporates the Hiller rotor head system (Figure 
2.1). The rotational inertia of the UAV main rotor is too much for the small servo 




Figure 2 . 7 Hiller Rotor System 
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actuators to overcome when an aircraft attitude change is demanded. Therefore, the cyclic 
control servos deflect the Hiller paddles, which either increase or decrease their angle of 
attack (AOA). As the AOA is changed, the lift on the paddles will increase causing the 
flybar to flap up or down (depending on the input). As the flybar flaps, it will change the 
orientation (tilt) of the swashplate causing a cyclic pitch change to the rotor blades. In a 
conventional helicopter, gyroscopic precession will cause any cyclic pitch input to be 
“felt” 90 degrees (or less) after it is input. For a perfectly “rigid” head or a teetering rotor 
system the precession lag is exactly 90 degrees following the input; for articulated rotor 
systems the precession is something less than 90 degrees dependent on the rotor flapping 
offset [Ref. 2], Similarly, the Hiller paddle experiences the same gyroscopic precession; 
therefore, the input is applied 90 degrees ahead of where we desire the control deflection 
to affect the rotor blade pitch. The blade azimuth position, X F, is zero over the tail and 
positive in the direction of rotor rotation [Ref. 2], In contrast to conventional U.S. 
helicopters that have a counter-clockwise rotating main rotor system, the main rotor on 
the Bergen Industrial rotates clockwise. In addition to the Hiller paddles, the rotor blades 
also experience a 90-degree precession; consequently, for a lateral or longitudinal cyclic 
attitude change, the actual input is applied 180 degrees ahead of the affected blade 
position [Ref. 3], For example, to input a nose down attitude change, the swashplate tilts 
forward which deflects the paddle at the = 90° position (negative AOA), the paddle 
flaps (down) to its minimum low position over the nose, at *F= 180°, decreasing the rotor 
pitch on the blade at the W = 90° position. As the blade continues to rotate to the W = 
180° position, it reaches its minimum flapping angle (flaps down), forcing the rotor disk to 
tilt forward causing a nose-down attitude. This nose-down attitude will result in an initial 
forward acceleration. 

C. EQUATIONS OF MOTION 

The behavior of the air vehicle in flight can be modeled as a combination of a 
number of interacting sub-systems (e g. the fuselage, main rotor, tail rotor, and 
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empennage). These sub-systems result in the aerodynamic forces and moments about the 
CG of the body (body axis system). Figure 2.2 shows the orthogonal body axis system 
used for aircraft dynamic analysis with the associated aerodynamic forces and moments. 




Figure 2.2 Helicopter Orthogonal Axis System [From Ref. 1] 

The equations governing these interactions are developed from the application of 
physical laws, i.e. Newton’s laws of motion relating the applied forces and moments to the 
resulting translational and rotational accelerations [Ref. 1], Where only the six DOF 
system is considered, the three translational velocities components are u, v, and w, and the 
three rotational velocities are p, q and r. The nonlinear equations of motion as presented in 
Reference 1 are given below. 

Force equations 

X 

ii = -{wq-vr) + — -g sin# (2.1) 

m 

Y 

v = -{ur - wp) + — + gcos0sin0 (2.2) 

m 

Z 

w = -(vp-uq) + — + qcos0cos</> (2.3) 

m 

Moment equations 

1 xx P ~ (Jyy ~ 1 zz I xz if PQ) 

I YyP ~ ZZ — I XX } r P I xz i r — P ) M 
I z zf = (I XX ~ I YY ")PQ I XZ iP — ^ 
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(2.4) 

(2.5) 

( 2 . 6 ) 



where the external forces (X, Y, and Z) and moments ( L , M, and N) are written as 
the sum of the contributions from the different aircraft components. Ixx, In, etc. are the 
fuselage mass moments about the CG, and m is the aircraft mass. The Euler angles, (j>, 0 
and % represent the orientation of the fuselage body axis with respect to an earth-fixed 
coordinate system. The external forces and moments, and Euler angles (save the yaw 
angle, I F) can be written in first order vector form: 

^ = f(x,u,t) (2.7) 

where x(t) is the column vector of state variables, x={u, w, q, 0 v, p, <f>, r } and 
u(t) is the vector of control variables. The control vector has four components: 
longitudinal cyclic, collective, lateral cyclic and pedals (directional control), 
U={5e,5c,5 a ,5p}. 

D. LINEARIZATION 

Using small perturbation theory, we assume that the helicopter behavior can be 
described as a perturbation from the trim condition. Assuming that the external forces and 
moments can be represented as analytical functions, we can expand the functions into a 
Taylor series about an operating point, and retain only the linear terms. The linearized 
equations of motion for the full six DOF system, describing perturbed motion about a 
general trim point, can be written as [Ref. 1]: 

x’ = Ax+Bu (2.8) 

where A and B are the system and control matrices respectfully, derived from the 
partial derivatives of the nonlinear function,/ (alternately referred to as F and G), where 
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and 




(2.9) 




( 2 . 10 ) 



The elements of the system and control matrices can be found in Appendix A. 
Note: the linear representation is valid only if the initial angular velocities are zero. 

E. STABILITY DERIVITIVES 

The elements of the A and B matrices are known as the stability and control 
derivatives. There are 36 stability derivatives and 24 control derivatives in the standard six 
DOF system. These derivatives represent the slope of the forces and moments at the trim 
point reflecting the strict definition of the stability and control derivatives [Ref. 1], The 
derivatives are non-dimensionalized by dividing the change in forces by the mass of the 
aircraft, and dividing the change in moments about each axis by the appropriate moment 
of inertia. A complete listing and explanation of all derivatives is presented in References 1 
and 2. 

F. SOFTWARE (JANRAD) 

System identification is the process of constructing a simulation model and 
associated parameters from experimental data. Numerous system identification programs 
are available to assist in the UAV “plant” identification process. Many of these programs 
use a curve fitting process by which experimental data is used to determine the system 
stability derivatives. To obtain the experimental data necessary to perform these functions 
requires a fully instrumented air vehicle. The lack of this capability, time requirements, and 
monetary constraints, led to the use of the JANRAD program to accomplish the task of 
system identification for the UAV. 
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The Joint Army/Navy Rotorcraft Analysis and Design (JANRAD) program was 
developed at NPS as an interactive, MATLAB' based program to meet the needs of the 
aerospace engineering student in preliminary helicopter design [Ref. 4], The JANRAD 
program is capable of performing both performance and stability analysis for any 
rotorcraft. The Performance program, described in the reference, determines the trim 
solution, and various performance parameters at a given flight condition necessary in the 
preliminary design. The Stability program calculates the stability derivatives for a given 
flight condition, and determines the state-space linear model at any trimmable point 
[Ref. 5], The stability derivatives are determined by using closed-form solutions whenever 
possible, or by solving multiple trim solutions about a nominal position. For the hover 
analysis, the stability derivatives from the main and tail rotors have the largest influence on 
aircraft response. 

For the short-term anticipated application requirements of the simulation model, 
the functional fidelity obtained from the JANRAD linearized model was deemed 
satisfactory. The JANRAD program was developed for full-scale helicopter design; 
therefore, the software code required slight modifications to accommodate the parameters 
representative of a scale-size helicopter model. Output fields were expanded to indicate 
values representative of the small UAV. Additionally, the output files and plotting routines 
were changed to allow compatibility with MATLAB upgrades. Continued software 
enhancements will make the JANRAD program more user-friendly, and increase the 
software’s analysis capabilities [Ref. 6], 

G. SIMULATION MODEL 

From the study of modem control theory, the NPS student has a working 
knowledge of the SIMULINK® simulation software. Therefore, its use was a natural 
progression when selecting simulation software. Other programs, such as the Systems 
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Build (XMATH) programs are available, but require an extensive time investment to 
become proficient. 

The basic SIMULINK Model developed uses state-space analysis of the system 
discussed above (Equation 2.8). The model in presented in the figure below (Figure 2.3). 




Figure 2.3 Simulink Block Diagram Model 

The analysis uses the F & G matrices input about a known trimmable condition. 
The F&G matrices were computed using the JANRAD Stability program. The matrix 
parameters were entered into a MATLAB' file, where the Simulink program was able to 
read the data. The input vector was designed such that single-channel inputs are possible 
while holding the remaining controls in the trimmed position. Recall: the control inputs are 
only changes (i.e. 5 e , 8 C , 5 a or 6 P ) from the trimmed condition. The outputs from each 
element of the state vector are presented with time history plots, such that any output 
from the state can be analyzed. From the output plots, dynamic response to control input 
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can be observed in the axis of interest (input), while also observing the cross-coupling 
effects experienced by the UA V. Information of UA V yaw angle, *F, is available by simple 
integration of the yaw rate, r. Additionally, acceleration data can be observed and later 
compared with flight test acceleration outputs. Prior to integrating the state-vector to 
obtain the axial velocities, u, v, and w, the state accelerations, w,v, and w, were 
obtained. 
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m. INPUT PARAMETERS 



The parameters required to conduct the performance, and stability and control 
analysis using the JANRAD program were determined either from measurement, table 
look-up or experimental tests. References (4) and (5) present a full listing of the required 
input variables. Reference (7) presents scaling as an alternative method to obtain aircraft 
physical parameters available only by experiment or flight data. The scaling method yields 
approximate results useful for the initial model development until more accurate means 
can be employed. A complete listing of all input data is presented in Appendix B. 

A. PHYSICAL MEASUREMENTS 

Obviously, the most desirable method is that of direct physical measurement. 
Measurements of those input variables that could be measured from the UAV were taken 
from an arbitrary reference datum. The datum, located roughly at the unmodified UAV 
CG, was designated as waterline (WL) zero, and buttline (BL) zero. The nose of the 
helicopter was designated as fuselage station (FS) zero. Additional input data that required 
“alternative” methods of collection are listed in the following sections. 

B. MOMENTS OF INERTIA 

When developing a mathematical model for an air vehicle, one cannot overlook the 
contributions of the aircraft’s mass-moments of inertia to its controllability and dynamic 
responses. The mass-moments of inertia represent the vehicle’s resistance to acceleration 
or rotation given a control input or external perturbation. Direct calculations of the 
moments of inertia can be accomplished by multiplying the mass of each component by the 
square of the distance to the body axis of rotation. For the scale UAV, this method is 
impractical because the air vehicle’s individual parts are too small and light to yield 
anywhere near accurate results [Ref. 7], Therefore, the aircraft’s moments of inertia must 
be determined by experimental methods. The reference describes two methods of 
developing a compound pendulum for obtaining the necessary moments experimentally: 
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the overhead pivot and the pivot at CG. The overhead pivot method of the compound 
pendulum was utilized. 

By suspending the UAV with small wires to a single pivot point on the ceiling, we 
can develop a compound pendulum (system) Knowing the vehicle mass, m , and distance, 
d, from the UAV CG to the pivot point, we can determine the moment of inertia about the 
helicopter’s CG, Icg, using the parallel axis relation: 



I = I CG + md 2 



(3.1) 



By giving the UAV a gentle push in a particular direction (along a body axis), we 
can oscillate the system, exciting the rotation of the body. The oscillatory period is 
determined by counting the number of cycles for a particular elapsed time. A cycle is 
defined as one complete oscillation to and fro. The period of the system will be dependent 
upon the body’s moment, Icg, the distance to the point of rotation, d , and the mass, m. 
The oscillatory period, Pm+s, is simply the total elapsed time (in seconds) divided by the 
total number of cycles. Reference 7 provides the equation for calculating Icg of the 
oscillating model: 



I 



CG 






-to 



(3.2) 



where W is the weight of the model, / is the distance from the pivot to the body C.G., and 
/o is the moment contribution of the supporting structure. To obtain the desired moments 
about all three body axes, the model was hung three different ways in order to obtain 
rotation about the axis of interest. Initially the tests were conducted using chains to mount 
the model to the ceiling. The results were suspect most likely due to the weight of the 
chain dominating the calculations leading to results, -1< Icg <1 [slug-ft 2 ]. Therefore, to 
eliminate the influence of the supporting structure, the experiment was repeated using 
lightweight, monofilament fishing line (/o assumed zero when using light-weight wire 



16 



supports). The results of Trial 2 show that the moment about the x-axis and y-axis are the 
same Since there was limited rotational motion about the axis of interest by displacing the 
system only a small amount, the system reacted similar to a point mass confirming these 
results. 

The experiment was conducted yet again (Trial 3) with a shorter distance from the 
UAV CG to the pivot point. This was done to further excite the rotational motion about 
the axis of concern. The results from this trial more closely match the values calculated 
from the scaling method (with a scaling factor, X=10.6). A complete listing of the 
experimentally determined mass-moments of inertia are presented below (Table 3.1). The 
experimental specifications and resulting calculations are presented in Appendix C. 



Trial 


Ixx [slug-ft 2 ] 


Iyy [slug-ft 2 ] 


Izz [slug-ft 2 ] 


Scaling 


.0374 


.2989 


.2615 


1 


.876 


1.32 


(.93) 


2 


.658 


.658 


1.887 


3 


.0709 


.3967 





Table 3 . 7 Moment of_ Inertia Results 



C. ROTOR FLAPPING MOMENT 

The rotor flapping moment is the mass-moment of inertia of the blade about the 
flapping hinge. The rotor flapping moment influences the rotor blade’s ability to flap due 
to blade pitch changes caused by cyclic inputs. The moment of inertia is defined as 
[Ref. 2]: 

R 

I = ^mr 2 dr [slug-ft 2 ] (3.3) 

0 
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where m is the specific mass of the blade (slug/ft), r is the radius of the blade element, and 
R is the total blade radius. For the main and tail rotors, the blade mass distribution was 
assumed uniform, therefore, the above relation can be simplified: 



where in this relation, m is the total blade mass (slugs). The results of these calcualtions 
are shown below. 



D. EFFECTIVE HINGE OFFSET 

To eliminate the rolling moment characteristic of the early autogyro and rotorcraft 
using rigid blades, the blade was allowed to flap about a flapping hinge. For rigid rotor 
systems, this flapping is accomplished by substituting a flexible section next to the hub 
[Ref. 2], The hinge offset of a fully articulated rotor system is merely the radial distance 
from the rotor hub to the rotor-flapping hinge. For a rigid rotor system, such as that used 
on the Bergen UAV, an effective hinge offset must be determined. The effective hinge 
offset is dependent upon the main rotor blade’s natural frequency, and rotor rotational 
velocity. The natural frequency of the blade is dependent upon the blade stiffness, blade 
length and mass. The stiffness coefficient of the blade is not a parameter normally 
specified in a RC helicopter owner’s manual; therefore, it must be determined 
experimentally. The basis of the experimental determination of the stiffness coefficient, 
El, of the blade is that the deflection, vmax, of a cantilever beam subjected to a point load, 
P, as illustrated in Figure 3 . 1 , is inversely proportional to El. 

This relationship is provided in Reference 8 by the equation: 




(3.4) 



lMRf=.034 [slug-ff 2 ] 



IxRf- 2.8 x 10‘ 5 [slug-ff 2 ] 




(3.5) 
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Figure 3. j Deflection of_ Cantilever Beam to a Point Load [From Ref. 8] 



where L is the beam length. To simplify the calculations, the blade was assumed 
homogeneous, and of constant cross section. The blade was placed in a vice, and using a 
fish scale, various loads were applied at the tip of the blade, and the deflections were 
recorded. Additionally a second method was performed where the tip was deflected a 
particular distance, and the applied load measured. 

The experimental data is plotted below (Figure 3.2) with calculations and results 
presented in Appendix D. With the blade stiffness coefficient determined, we can now 
calculate the natural frequency of the blade, <**,. Reference 9 provides the required relation: 



where co \ is the natural frequency for the first bending mode of the blade, and m is the mass 
of the blade. 




(3.6) 
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Figure 3.2 Plot of Beam Deflection vs. Applied. Load 

From the blade natural frequency and the rotational frequency of the main rotor, Q, we 
can determine the effective hinge offset ratio from the equation provided by Reference 2 
(pg. 457): 




(3.7) 



From this analysis, the frequency ratio, co n /Q <0 leads to a value of (e/R) e ff < 0. 
The physical implications would be that the rotor system exhibits no effective hinge offset. 
As was described in Chapter 2, the Hiller paddles drive the cyclic pitch of the rotor blades, 
and are free to flap about the rotor head. The Hiller paddles essentially operate the same 
as a teetering rotor system which has a hinge offset of zero [Ref. 2]. With a hinge offset of 
zero, the rotor flapping response is exactly 90° out of phase with the applied cyclic pitch. 
Which agrees with the basic design of the UAV cyclic control phase relationship where the 
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swashplate input is exactly 90° out of phase with the desired disk movement (from the 
action of the Hiller flybar). 

E. AERODYNAMIC INPUTS 

Several of the required input data for the JANRAD program were neither available 
by direct measurement nor experimentally; therefore, several resources were required to 
satisfy the required data input fields. Many of these data were aerodynamic coefficients of 
control surfaces, or pressure ratios experienced by these aerodynamic surfaces used in the 
performance analysis. Below is a listing of these input data with a brief description, 
reasoning and reference source. 

Flight Conditions . Required input data concerning the flight conditions were 
input as a standard flight test day of 60 degrees Fahrenheit and pressure altitude of 100 ft 
MSL. For comparison with actual flight data, flight condition information is available from 
the local airport weather center. 

Forward Velocity : For the hover analysis, forward velocity is equal to zero 
(no wind condition). If performing an analysis of forward flight, UAV airspeed data is 
required in order to conduct an accurate simulation analysis. At present, this capability 
does not exist from the vehicle, however, a simple method can be devised by determining 
the time to fly between two points a known distance apart yielding an approximate, no 
wind, velocity. Future sensor enhancements may incorporate an internal capability. Note: 
the JANRAD analysis does not recognize forward airspeeds below 1 2 knots, and accuracy 
is limited below 50 knots [Ref. 4], 

Fuselage Downwash Ratio (due to rotor) : The fuselage downwash ratio 
corrects for the interference of the fuselage due to the downwash of the main rotor 
system. Using the analysis in Reference 2 (Figure 8.11 pg. 494) with XVR= -1.07 and 
ZVR= +0.19, an approximate value of vh/vi of 1.5 was selected. 

Main and Tail Rotor Lift Curve Slope . Exact airfoil data for the Bergen 
UAV was unavailable from the manufacturer, and time constraints limited windtunnel 
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testing. References 2 and 10 provide example lift curve slope values for helicopter rotor 
blades. From the values presented, a typical value of 5.73 (per radian) was used for the 
analysis. 

Rotor Rotational Velocity; . The rotational velocity of the main rotor system 
was determined using a RPM checker. With the UAV hovering at its operational rotation 
velocity, the rotational velocity was collected in rotations per minute (RPM). These values 
were verified during the flight tests by analysis of the vibration data. The rotational 
velocity of the tail rotor was determined by using the gearing ratio between the main rotor 
and the tail rotor (1: 4.6). The throttle controller on the UAV is not a constant RPM 
system, i.e. as rotor loads increase, the rotor rotational velocity increases. Limited data 
was collected during initial flight-testing; therefore, the operating RPM should be verified 
during each flight prior to performing the computer simulation. 

Rotor Blade Airfoil '. As stated above, the UAV airfoil data was 
unavailable. An additional limitation on the analysis performed is that the JANRAD 
software used contains only data for three airfoil models: NACA 0012, Boeing VR-12 and 
the HH-02 [Ref. 5], The main rotor has a cambered airfoil; therefore, the analysis was 
performed using a cambered airfoil. The VR-12 was selected. 

Honzgntgl and Vertical Tail Coefficients of Lift and_ Drag : Both the 
horizontal and vertical tails are flat plates. Airfoil data for these surfaces is available in 
References 2, 11 and 12. For a flat plate airfoil, a value of the maximum lift coefficient, 
Cimax, equal to 0.8 can be expected [Ref. 2], The coefficient of drag, O, for a flat plate is 
between 0.004 and 0.006 [Ref. 12]; a value of 0.005 was used. 

Horizontal and Vertical Tail Lift Curve Slope . The expected values of the 
empennage surface lift curve slope are dependent upon the airfoil sections and their aspect 
ratios [Ref. 2 ]. For the flat plate, representative values are plotted against effective aspect 
ratio, ARe. Using the AR* of the horizontal surface, the value of dC\/da was obtained. 

Dynamic Pressure Ratio . The dynamic pressure ratio, qn/t 7, provides for a 
corrected flight condition (dynamic pressure) experienced by the empennage due to the 
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presence of the main rotor system (rotor wake). From the curves provided in Reference 2, 
an approximate value of 0.6 was obtained. 

Rotor Downwash Ratio : The rotor downwash ratio, v H /v t , is a ratio of 
vertical velocities where v H is the vertical velocity at the horizontal stabilizer, and vj is the 
vertical velocity of the induced velocity in the plane of the rotor. This ratio is dependent 
upon the horizontal and vertical position of the stabilizer with respect to the rotor hub. 
From the data curves provided in Reference 2, a representative value of 1.5 was used. 

Fuselage Downwash Ratio : The fuselage downwash ratio gives the 
coefficient of the downwash effect of the fuselage on the horizontal stabilizer. Reference 2 
provides coefficient values for use as an input parameter for helicopters with or without 
wings. The original UAV configuration has no wing; therefore, the value of dSF/dctF 
selected was 0.06. 

F. CONTROL RIGGING DIAGRAMS 

To determine the required data for computing the control derivatives, the UAV 
control-rigging scheme was needed. For the sample helicopter used in Reference 2, the 
rigging charts are plotted with degrees of rotor blade pitch, either A\ or B\ for cyclic 
inputs, versus inches of cyclic stick deflection. Because the UAV is controlled remotely 
using radio inputs, or pulse width modulation (PWM), it was desired to plot control 
surface movement versus PWM. Using this methodology, we can determine the applied 
control input by comparing the “captured” PWM signal sent to the UAV during a certain 
maneuver, then model the input and response for comparison with flight data. 

The rigging data was collected separately from each of the four inputs: longitudinal 
cyclic, collective (throttle), lateral cyclic, and directional pedals. For each measurement 
taken, the helicopter and flybar were leveled both longitudinally and laterally using a 
bubble-sight level. Using the Schluter AOA Gauge (Figure 3.3), we were able to measure 
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Figure 3.3 Schluter Blade Angle Measurement Device 

rotor blade AOA while the corresponding PWM signal required to maintain that control 
position was collected (at the ground station). For longitudinal cyclic, and collective pitch, 
blade AOA was measured with the blade at { F= 90° (90° ahead of desired reaction due to 
gyroscopic precession). Lateral cyclic data was collected with the rotor at the y F= 180° 
position. When taking AOA measurements in both cyclic channels, it was necessary to 
apply PWM collective control to bring the neutral cyclic pitch to 0° AOA. Figures 3.4 
through 3.7 show the experimentally collected data. A trendline was used to determine the 
slope of the data for input to the control matrix and compensate for measurement errors. 
During the simulation run, the control inputs applied will be with respect to PWM signal 
rather than inches of control stick deflection. 
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Longitudinal Cyclic Rigging Curve 



4S00- 




, , — 1200 - 1 t . 

- 8-3 2 7 

Cyclic Pitch, B^deg) 



Figure 3.4 Longitudinal Cyclic Rigging. Curve 
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Figure 3.5 Collective Rigging Curve 
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Lateral Cyclic Rigging Curve 




Figure 3. 6 Lateral Cyclic Rigging. Curve 



Directional Control Rigging 




Figure 3 . 7 Directional Control Rigging Curve 
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IV. SIMULATION ANALYSIS 



A. JANRAD IN HOVER 

The first step in the analysis of our model is the use of the JANRAD Performance 
program to evaluate the UAV trim condition. Once the helicopter trim position is 
determined, the trim solution can then be used to evaluate certain performance parameters. 
The JANRAD program was developed as a preliminary design tool, and uses the harmonic 
balance method to trim the rotor using actual 2D airfoil data. The procedures for running 
the original JANRAD program are outlined in Reference 4. The enhanced JANRAD 98 
version is also available with increased capabilities and a user-friendly Graphical User 
Interface (GUI) [Ref. 6], 



1. Performance Calculations 

The analysis was conducted using the input data developed in Chapter 3 with a 
complete listing in Appendix B. The results of the Performance analysis are listed below 
(Table 4.1) in the output format from JANRAD. Many of the values are zero with the 
UAV in a hover. 



Fuselage drag = 0.00 lbs. 

Rotor drag = 0.00 lbs. 

Wing lift = 0 lbs. 

Wing drag = 0 lbs. 

Horizontal tail lift = 0.00 lbs. 
Horizontal tail drag = 0.00 lbs. 
Vertical tail side force = 0.00 lbs. 
Vertical tail drag = 0.00 lbs. 

Tip path angle = 0.00 degs 
Rotor coning angle = 1.58 degs 
Location of mean thrust (r/R) = 0.75 
Collective pitch at .7 r/R = 4.61 degs 
1st lat cyclic term-Al (deg) = 0.00 



1st long cyclic term-Bl (deg) = 0.00 

Solidity (sigma) = 0.045 

Disk loading = 0.78 lbs/ft A 2 

Figure of Merit = 0.46 

CT/sigma = 0.034 

CQ/sigma = 0.0021 

CH/sigma = 0.0000 

Tip mach of the adv blade = 0 415 

Advance ratio = 0.000 

Rotor thrust required (TPP) = 19.67 lbs. 

Rotor power required = 1 .01 h.p. 

Rotor torque = 3.37 ft-lbs. 



Table 4. j Listing, of JANRAD Performance Output Data 
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Verification of the JANRAD calculations is only possible by generic computations 
from Reference 2. Testing equipment is also not available for comparison or verification. 
The trim solution determined in the Performance program is also used as a starting point 
for the Stability Analysis program. 

2. Stability Analysis 

The JANRAD Stability Analysis Program [Ref. 5] provides the stability and 
control derivatives (F and G matrices) which become the plant of our simulation model. 
The calculated matrices are listed in Appendix E. These matrices can also be used to 
perform an open-loop stability analysis by analyzing the open loop eigenvalues, and 
performing a frequency response analysis or “Bode Diagram”. 

For small amplitude stability analysis, helicopter motion can be considered to 
comprise a linear combination of natural modes, each having its own unique frequency, 
damping and distribution of the response variables [Ref. 1]. Using the stability matrix, F 
(or A), we can describe the free motion of the UAV in the form: 

x - Ax - 0 (4.1) 

subject to the initial conditions, x(0)=xo. The natural modes are described as 
linearly independent so that no single mode can be made up of a linear combination of the 
others. Then, if a single mode is excited, we should expect the motion to remain in that 
mode only. In the table below (Table 4.2), we provide a listing of the system eigenvalues, 
which are also plotted in Figure 4.1. From the complete list of eigenvalues, we can 
decouple the system in order to separate the eigenvalues into two sets, the longitudinal 
and the lateral modes. The characteristic roots, which appear as complex conjugate pairs, 
represent an oscillatory response for that mode. 
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Eigenvalues 


Freq. (rad/sec) 


Damping 


-5.2075 


5.2075 


1.00 


-1.8853 


1.8853 


1.00 


- 0.7560 


0.7560 


1.00 


-0.6625 


0.6625 


1.00 


0.2149+0.581 li 


0.6196 


-0.3469 


0.2149-0.581 li 


0.6196 


-0.3469 


0.1890+0.4094i 


0.4510 


-0.4191 


0.1890-0.4094i 


0.4510 


-0.4191 



Table 4.2 Eigenvalues o f the Linearized Model 



Root Locus 




Figure 4 . 7 Open-Loop Eigenvalue Plot 
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The stability of the UAV can be discussed in terms of the individual eigenvalues, 
which is determined by the sign of the real part of the root. A root with a negative real 
part, or located on the left-hand side of the argon plane is stable, while a root with a 
positive real part (right-hand side of plane) is unstable. 

We can see there exists two unstable, oscillatory roots. The first pair: 

X=0.1890 ±0.4094i 

correspond to the longitudinal phugoid, or long period, mode. This mode also has a 
negative damping value which can lead to a divergent longitudinal mode when excited. 
The long period (P=13.93 seconds) characteristic of this mode can be easily controlled by 
a pilot in the loop, or automatic flight control system. The second pair of unstable, 
oscillatory roots: 

A.=0.2149±0.581 li 

correspond to the lateral dutch roll mode. This mode has a period of 10.14 seconds, and 
again could be pilot or stability augmentation system controlled. For the open loop 
simulation analysis, these instabilities may prove problematic. 

An additional method for analyzing the relative stability of a system is by the use of 
frequency-response tests. By the term frequency-response, we are referring to the steady- 
state response of the system to a sinusoidal input. Using this method, we vary the 
frequency of the input to the system over a certain range and study the resulting system 
response. A complete frequency-response analysis would include the response of each 
state of the system to all the possible inputs. For the helicopter, this would involve 
analyzing the eight states of the system to each of the four possible inputs. At this point, 
we are not concerned with the effects of cross-coupling (off-axis response to control 
input); therefore, we will concern ourselves with looking at the responses along the three 
orthogonal axes from inputs along those same axes (i.e. vertical response to a vertical, or 
collective, input, etc.). The frequency-response analysis is graphically depicted by either 
the Nyquist Plot or the Bode Diagram. To plot system response magnitudes, the Bode 
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diagram uses a semi-logarithmic plot of system response in decibels, dB, plotted against 
input frequency in radians. A second plot shows the system phase angle — the angular 
difference between the harmonic input and the system response. One important 
characteristic we can determine from the Bode plot is the “bandwidth”. The bandwidth is 
indicated by the frequency at which the response has decreased by 3 dB, or about 71 
percent of the initial response. 

By applying any control input, we can determine the frequency-response to the 
mathematical model by using the MATLAB command “bode (F, G, C, D, iu)”, where F, 
G, C and D are the matrices of the state-space model, and iu produces an input to the zth 
element of the control input vector. Figure 4.2 provides the Bode plot for a longitudinal 
response (forward velocity, u) to a longitudinal cyclic input, 5 e . The effective bandwidth 
for the longitudinal cyclic is approximately 0.50 rad/sec. From the lower plot, we note that 
the system displays a phase “lead”. This is characteristic of systems with an unstable root 
in a particular channel. Figure 4.3 is the result of applying a collective input to the model, 
and analyzing the frequency-response in the vertical axis. The bandwidth for this input is 
approximately 1.5 rad/sec with a phase “lag” as indicated by the lower plot. Applying a 
lateral input to the system (Figure 4.4), and looking at the frequency-response in the 
lateral axis yields similar results as seen in the first plot (Figure 4.2). The gain drops 3 dB 
at a frequency of approximately 0.70 rad/sec. The phase shift again is shown to be 
unstable. 
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Longitudinal Cyclic to Forward Velocity 





Figure 4.2 Bode Plot of Longitudinal Cyclic to Forward Velocity, 
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Phase deg Gain dB 
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Collective to Vertical Velocity 




Figure 4, 3 Bode Plot oj_ Collective to Vertical Velocity. 
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Lateral Cyclic to Sideward Velocity 





Figure 4. 4 Bode Plot of Lateral Cyclic to Sideward Velocity 
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B. SEMULINK MODEL IN HOVER 



The above analysis provides some insight to the inherent stability of the simulation 
model by using the frequency-domain method. By applying discrete inputs to our 
simulation model, we can analyze the responses in the time domain. 

Using our Simulink model developed in Chapter 3, we can input either a unit step, 
impulse or doublet to perturb the model from its trim condition, and examine any or all of 
the eight “states” of the system. As was done in the previous section, our initial analysis of 
the system will be primarily restricted to analyzing the response along the axis of the 
applied input. Recall when we plotted control rigging schedules for the four UAV inputs, 
each was plotted versus the radio controller signal, or PWM, required to maintain the 
control position. Therefore, to apply a realistic input to the simulation model, the 
magnitude of the input signal was adjusted to change the control surface by approximately 
one degree (either cyclically or collectively). Small inputs were used in an effort to predict 
the actual UAV responses (displacements) anticipated during testing to remain within the 
flight-test operating area. The following time-history simulations show both short-term 
initial UAV responses, as well as “long-term” oscillations. 

Figure 4.5 shows the applied longitudinal doublet input and the resulting forward 
acceleration, velocity, and pitch angle response. Recall from the stability analysis that the 
longitudinal phugoid mode was unstable, causing extreme divergent oscillations as time 
progresses. Of particular interest is the period of the phugoid oscillations, from the figure, 
we can see an oscillation period of approximately 14 seconds which matches our earlier 
analysis of the characteristic roots. Because of the large scale along the acceleration 
vertical axis, it is difficult to discern the magnitude of the initial acceleration response; 
therefore, the simulation was replotted showing the initial, short-term response (Figure 
4.6). The initial longitudinal acceleration shows an acceleration of only 0. 1 G’s. 
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Simulation Longitudinal Doublet 




Figure 4.5 Long = Term Longitudinal Doublet Simulation 
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Figure £ 6 Short-Term Longitudingl Doublet Simulation 
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A similar input was applied to the vertical axis by applying a collective doublet to 
the simulation model, and plotting the acceleration and yaw angle, *F. Figure 4.7 indicates 
the model responded as expected; the model accelerates up and down following the 
collective inputs, then, because the system is stable in the vertical channel, returns to the 
trim position. The model yaws in response to the collective input, and also returns to the 
trim condition. At approximately 10 seconds in the simulation, the yaw angle starts to 
deviate once again in response to coupling instabilities in the other axes. 

For the lateral mode, a lateral (cyclic) impulse was applied; the impulse and 
doublet are both useful in exciting the long-period instabilities of the system. The resulting 
simulation is shown in Figure 4.8. Similar to Figure 4.5 above, the system eventually goes 
divergent due to the unstable dutch-roll characteristic root. The period of this unstable 
root from above was calculated to be approximately 10 seconds, which is depicted in the 
resulting acceleration, velocity and roll angle shown in the figure. The short-term response 
(Figure 4.9) indicates an initial lateral acceleration of approximately 0.2 G’s. 
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Figure 4 . 7 Collective Doublet Simulation 
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Figure. 4.8 Long = Tenn Lateral Impulse Simulation 
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Figure 4.9 Short-Term Laterdj Impulse Simulation 
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V. FLIGHT TESTING 



A. INSTRUMENTATION OF FLIGHT VEHICLE 

To obtain useful flight data to validate our simulation model, the UAV requires the 
incorporation of onboard instrumentation. Ideally, access to UAV airspeed, altitude, 
attitude, and acceleration would be available to accurately reconstruct air vehicle 
responses experienced during the test flight. For the initial flight-testing phase of the 
research, the goal was to obtain acceleration data in an out-of-ground effect (OGE) hover. 
The original instrumented configuration had an inertial measuring unit (IMU) with 
associated downlink telemetry devices mounted under the UAV on a modified landing 
gear with avionics rack. Reference 13 provides a detailed description of the avionics suite 
design. 

The IMU, also referred to as an attitude, heading reference system (AHRS), 
provides bank (roll angle), elevation (pitch angle), heading (yaw angle), translational 
accelerations, and rotational accelerations along the X, Y and Z orthogonal axes. Due to 
space constraints, it is usually impossible to place the IMU at the ideal location, the air 
vehicle CG; therefore, the output data from the EMU must be transformed from its 
position to the aircraft’s CG to the inertial reference system. Other onboard sensors can 
also be incorporated into the avionics suite to improve the collection of UAV response 
data, such as airspeed indicators and global positioning system (GPS) sensors. 

B. PHYSICAL PARAMETER MODIFICATION 

With the addition of instrumentation to the UAV, the performance and dynamic 
responses of the air vehicle will be changed due to the increased weight (mass) and altered 
mass distribution. Added weight will not only shift the CG position, but will change the 
mass-moments of inertia. Since the simulation model was developed about a ‘trim 
condition’, we must repeat the stability analysis to reflect the changes to the system’s 
dynamic characteristics. In order to recalculate the new trim position, we must start by 
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changing the input parameters of the JANRAD program to obtain a new set of system and 
control matrices (F and G). In the following sections, we explain some of the methods 
used in determining the ‘new’ input variables. 

1. Center of Gravity 

The addition of onboard sensors and telemetry gathering equipment will alter the 
weight and mass distribution of the air vehicle. It is therefore necessary to make 
accommodations to our simulation model to reflect these changes experienced by the 
UAV. In addition to changing the vehicle weight, we must account for the altered mass 
distribution. The position of each additional component was determined with respect to 
the original UAV fuselage station (FS), waterline (WL) and buttline (BL). The component 
mass was determined by using a calibrated scale and dividing by the gravity constant. An 
Excel' spreadsheet was developed to automatically account for the change in weight and 
CG position (Appendix F). The spreadsheet can easily be updated by simply annotating the 
addition of the added component on the worksheet; the component mass and location data 
are already determined. A second work sheet, discussed below, computes the change in 
UAV mass moment of inertia as a result of the new mass distribution. 

2. Moments of Inertia 

The change in mass distribution will alter the dynamic response of the UAV when 
perturbed from the trim position. A good example of the change in the moment of inertia 
is the spinning ice skater; with the hands extended, the skater spins slow, but when the 
hands (mass) are brought in close to the body (axis of rotation) the moment of inertia is 
reduced, and the rate of spin increases. Therefore, to accurately predict the dynamic 
response of the UAV to a control input, the moments of inertia will need to be 
recomputed accounting for the added mass and mass distribution. Using the Parallel Axis 
Theorem, the Excel program described above will automatically compute a new moment 
of inertia about the updated center of gravity (page 2 of Appendix F). For simplification, 
each added mass was treated as a point mass, i.e. the moment of inertia about its own CG 
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was neglected. The standard position of each component was located on the avionics rack 
with respect to the rack forward centerline. 

After the moment of each component was determined, the new moment 
about the desired axis was obtained by using the relationship that the total moment of 
inertia of a body is equal to the sum of the individual parts about a common reference, i.e. 
the updated center of gravity (CG). 

3. JANRAD Input Modifications 

To obtain an updated set of ‘plant’ and control matrices from the JANRAD 
program, the process described in Chapter 3 will be repeated using the data obtained for 
the new UAV configuration. In addition to changing the CG position and mass-moments, 
additional parameters will need to be updated prior to running the analysis. Table 5.1 
below provides a listing of all JANRAD input variables to be modified. The values used in 
the analysis are listed in Appendix B. 



Parameter 


Units 


JANRAD Input Screen 


Forward Velocity 


[knots] 


JANRAD “Edit Menu” 


Temperature 


[°F] 


JANRAD “Edit Menu” 


Pressure Altitude 


[ft] MSL 


JANRAD “Edit Menu” 


Equivalent Flatplate Area 


rft 2 i 


JANRAD “Edit Menu” 


Vertical Projected Area 


[ft 2 ] 


JANRAD “Edit Menu” 


Gross Weight 


Pounds [lbs ] 


JANRAD “Edit Menu” 


CG Position 


[ft] 


Stability & Cntl (3 of 3) 


Mass-moments of Inertia 


[slug-ft 2 ] 


Stability & Cntl (3 of 3) 


Main Rotor Rotational 
Velocity, 


[rad/second] 


JANRAD “Edit Menu” 


Tail Rotor Rotational 
Velocity, Qtr 


[rad/second] 


Stability & Cntl (1 of 3) 



Table 5 . 7 Altered Input Variables to JANRAD Program 
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C. REQUIRED MODIFICATIONS TO OUTPUT DATA 



The output telemetry collected during the flight-testing does not accurately 
describe the responses of the body unless the information is taken with respect to the body 
(UAV) CG. Therefore, we must convert this data to the CG position of the body, then to 
the inertial reference, for comparison with the simulation model. Additionally, dependent 
upon the noise attenuation capabilities of the acceleration equipment used for collecting 
data, a considerable amount of vibration noise can be present in the output data, masking 
the desired translational accelerations. 

1. Data Translations 

Using a fully instrumented air vehicle, we desire access to all the states of the 
system (i.e. velocities, rotational rates and attitude). To accurately reconstruct the 
responses of the helicopter, it is necessary to translate these responses from the position of 
the IMU to the CG of the body. The IMU is the desired type of measurement device for 
collecting this data because of the depth of output information available. Due to vibration 
problems experienced during the initial instrumented flights, the avionics rack with 
associated measurement devices was unserviceable for mounting on the UAV during the 
later test flights. As an alternate method of collecting acceleration data, a Crossbow 
three-axis accelerometer was temporarily mounted on the UAV to meet research 
completion requirements. It was necessary to tether, or hardwire, the UAV to the data 
collection ground station in order to collect the data from the sonar altimeter and 
accelerometer. With the vehicle tethered, the motion of the UAV was restricted to the 
limits of the umbilical. As mentioned above, we normally would use a transition matrix to 
transfer the data from the sensor to the CG of the air vehicle, but for this configuration, 
we must assume the accelerometer at the body CG. The translation of accelerometer data 
is not possible because the three-axis accelerometer does not provide angular rate 
information 
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2. Vibration Noise Elimination 



Unfortunately, the accelerometers utilized could not discriminate the desired 
translational accelerations from the high-frequency vibrations inherent of all helicopters. 
The major contributors of these additional vibrations are the main and tail rotors, and the 
helicopter power plant. The main rotor system acts like a “filter”, and only those 
vibrations with a frequency multiple of its rotational velocity (2P, 4P or 6P for a two 
bladed rotor) are translated through the system. Figure 5.1 (top figure) shows a sample of 
some raw acceleration data. It is difficult to identify the translational accelerations in 
response to pilot control inputs. Therefore, to make the acceleration data useful, it was 
necessary to filter out these extraneous, high frequency accelerations. Initially, a spectral 
density analysis was used to confirm the main rotor rotational frequency of 27 Hertz. With 
this information, the raw data could be filtered to eliminate the higher frequency noise 
from the desired acceleration data. Accelerometer data was filtered using a 4 th order 
Butterworth Low-pass filter with a cut-off frequency of 12 rad/sec. The filtered data 
(Figure 5.1 bottom plot) can now be used to identify the aircraft responses for analysis. 

D. FLIGHT DATA ANALYSIS 

The goal during the flight -testing phase was to replicate system responses available 
from the simulation model. The Simulink model described in Chapter 3 allows the system 
to be perturbed about a trim position by either a control step input, impulse or doublet 
(control input reversal). It was desired to have the UAV pilot replicate the identical type 
of inputs as were used in the simulations. This would simplify the comparison of UAV 
responses to those of the simulation model. Because of the limitations imposed by the 
“umbilical” used to transfer acceleration and altitude data to the ground station, the inputs 
to the UAV had to be minimized to maintain an acceptable margin of safety for the UAV. 
The test runs conducted included a series of longitudinal and lateral cyclic and collective 
doublets. Figure 5.2 (upper plot) shows the type of input used in the simulation trials. The 
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Figure 5 , l Comparison of Raw Acceleration Data with Filtered^ Data 
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Figure 5.2 Comparison of Simulation Input to Test-Flight Input 
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lower plot shows the actual PWM (corrected about a trim position) from an actual test 
flight It is apparent that comparing the responses of the model with the UAV using these 
different types of inputs would not be realistic. Therefore, in the comparison of the 
simulation model with the actual air vehicle required subjecting the simulation model to 
the same input as the UAV was subjected. The test flight PWM collected by the ground 
station was converted to a MATLAB compatible format, then input into the Simulink 
model for each simulation run. 

E. COMPARISON OF SIMULATION AND FLIGHT DATA 

The purpose of the flight-testing was the verification of the simulation model, and 
the validation of the Sonar Altimeter operation as the subject of Reference 13; therefore, 
the limited test time was divided between meeting both objectives providing limited test 
flight-data. Two test runs involved subjecting the UAV to collective doublets, while the 
lateral and longitudinal doublets were only examined on one trial each. Control input data 
and UAV accelerations were collected on separate data loggers which were not 
synchronized and used different sampling rates; therefore, the data had was synchronized 
by hand to yield an accurate comparison between simulation and test data. Many of the 
figures below are plotted versus sampling rate, or filtering rate vice time in seconds. All 
accelerations are given in “G’s” with altitude given in feet above ground level. 

The first mode we will analyze is the two vertical collective doublet test-flights. 
Figure 5.3 shows the actual collective PWM input from the test-flight and simulation, the 
simulation responses in G’s, the filtered acceleration data in the vertical axis and the 
altimeter data from the sonar altimeter. Figure 5.4 displays the same type of data for the 
second collective doublet flight. For both flights, it is apparent the simulation model tracks 
closely the actual responses of the UAV. Note: an increase in collective leading to a 
positive “G” and climb is indicated by a decrease in the PWM signal. 

The next series of flight-test data is for the lateral cyclic doublet. Figure 5.5 shows 
the lateral input, the simulation acceleration and the lateral test G’s loads. The simulation 
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Figure 5.3 Vertical Response to a Collective. Daublej (Flight 1} 
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Figure 5.4 Vertical Response to a Collective Doublej (Flight 5} 
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Figure 5.5 Lateral Response to a Lateral Cyclic Doublet (Longterm}. 



53 



model, because it is unstable in the lateral mode, yields divergent oscillations as time 
progresses. Therefore, the large scale used on the y-axis gives no indication of the initial 
lateral accelerations of the model. The period of the oscillations should be noted since it is 
on the order of 10 seconds, equal to the period of the unstable dutch-roll mode. To 
determine the initial response of the model due to the cyclic control input, we have 
replotted the same test data on a shortened time scale (Figure 5.6). From the figure, we 
can determine that the initial response of the simulation model is of the same acceleration 
magnitude as the test UAV. 

Subjecting the simulation model to a longitudinal cyclic doublet, the results were 
similar with those referred to above for the lateral mode. Figure 5.7 shows the long-term 
acceleration (G’S) and pitch angle (radians) response of the simulation model and UAV to 
a series of longitudinal doublets. Again we note the simulation model starts divergent 
oscillations due to the first longitudinal doublets leading to undetectable initial response 
data. Of note is the period of the acceleration and pitch angle oscillations; from Chapter 4, 
we identified the phugoid mode to have a period of approximately 14 seconds which is 
shown in the figure. The same flight data was plotted for a short duration (Figure 5.8) 
yielding initial responses of the simulation model similar to the initial responses, of the 
UAV. 
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Figure^ 5. 6 Lateral Response to a Lateral Cyclic Input (Short-Term) 
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Figure 5. 7 Longitudinal Response to a Longitudinal Cyclic Input CLonv-Term) 
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Figure 5.8 Longitudinal Response to a Longitudinal Cyclic Input (Short-Term) 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The rotary wing flight vehicle presents an interesting and complicated dynamic 
system for the stability and control analyst. To accurately model this complicated, non- 
linear, nine-degree of freedom (DoF) system requires complex algorithms and analysis. A 
common approximation of using a linear, 6 DoF system is useful for preliminary analysis 
using small perturbations about a trim condition. It was the effort of this research to use a 
linear analysis to develop a simulation model of the Bergen rotary wing UAV capable of 
predicting its responses to pilot controlled inputs. 

The linear model developed was capable of simulating the vertical accelerations 
experienced by the UAV along the vertical axis when subjected to the identical collective 
inputs from the flight-tests. Problems arose when the model was subjected to inputs along 
the lateral and longitudinal axes. The initial response of the model along these axes 
followed the UAV for response direction and magnitude of accelerations. But, over the 
long term, inherent instabilities of the UAV and simulation model led to divergent 
behavior when perturbed from the trim condition. The time period of these unstable 
oscillations was within the capabilities of the UAV pilot to control, but complicated 
matters for the simulation model. As these unstable perturbations persisted, additional 
inputs were made to maintain control of the UAV which further perturbed the model 
which was, at this point, no longer in a trim condition. 

The linear model could be used in the development of an altitude controller should 
the vehicle be used in a hover or in slow near-hovering flight. A more robust, non-linear 
model is needed to accurately predict the dynamic response of the UAV in other flight 
regimes. 
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B. RECOMMENDATIONS 



Because of the complexity of the rotary wing vehicle with its interacting 
components and cross-coupled dynamic responses, a non-linear model is presumably 
needed to accurately predict the UAV responses to controller inputs. To accomplish this, 
the UAV would require proper instrumentation providing access to all states of the 
system. With the collection of flight data, one of the many system identification programs 
could be used to identify the system, and modify the current model or develop a 
completely new one. 

To ease the complexity of the model responses, only a single test-flight control 
signal was input at any time. With a more robust model, it may be possible to apply all 
four inputs to the model from the control data collected during the test-flights. 
Additionally, an attempt was made to more closely model the actual UAV by the 
incorporation of a feedback stabilization yaw-damper. To facilitate remote-pilot operation 
of the UAV, a yaw rate gyro with feedback control is incorporated on the Bergen RC 
helicopter. This eliminated aircraft yaw during large power changes. With the feedback 
controller on the simulation model, the initial yaw angle experienced was reduced, but the 
inherent instabilities of the system were not eliminated. Applying feedback controllers to 
the unstable axes may increase the desired response of the simulation model. 

Although, all these capabilities currently exist, and are used to control a wide range 
of UAV platforms, this research, though not groundbreaking, provides the student a 
foundation for the development of simulation models and the control of flight vehicles. It 
is imperative to expand on the work presented in this report, and that presented in 
Reference 13 to fully realize the potential of past, present and future NPS endeavors in 
this area. 
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APPENDIX A: STABILITY AND CONTROL DERIVATIVES 
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Appendix A: Stability and Control Derivatives (A and B Matrices) 

Here is a complete listing of the stability and control derivatives that comprise the stability and control matrices (A and B, 
respectfully) [Ref 1], 



APPENDIX B. JANRAD INPUT DATA 



The following is a comprehensive listing of all input variables used in the JANRAD 
Performance and JANRAD Stability analysis programs. 



*** JANRAD Performance Input Data *** 



Forward velocity = 0 kts 

Temperature = 60 degs F 

Pressure altitude = 1 00 ft 

Gross weight = 19. 16 lbs 
Number of blades = 2 

Rotor radius = 2.83 ft 
Blade mean chord = 0.20 ft 
Blade twist = 0.00 degs 
Blade airfoil = VR-12 
Blade lift curve slope = 5.73 
Blade weight = 0.42 lbs 
Rotational velocity = 164.00 rads/sec 
Blade grip length = 0.27 ft 
Hinge offset = 0.00 ft 
Equivalent flat plate area = 0.67 ft 2 



Vertical projected area = 1.71 ft 2 
Wing area = 0.00 ft 2 
Wing span = 0.00 ft 
Wing CL = 0.00 
Wing CDo = 0.0000 
Wing efficiency factor = 0.00 
Horizontal tail area = 0. 12 ft 2 
Horizontal tail span = 0.66 ft 
Horizontal tail CL = 0.00 
Horizontal tail CDo = 0.0050 
Vertical tail area = 0.22 ft 2 
Vertical tail span = 0.92 ft 
Vertical tail CL = 0.00 
Vertical tail CDo = 0.0050 
Auxiliary thrust - 0 lbs 



*** JANRAD Stability Input Data *** 
*** INPUT DATA (screen 1 of 8) *** 



Flight Conditions 
Forward velocity = 0 kts 

Temperature = 60 degs F 

Pressure altitude = 1 00 ft 

Auxiliary thrust - 0 lbs 

Fuselage 

Gross weight = 19. 16 lbs 
Ixz = 0.000 slug ft 2 
Downwash ratio = 1.50 



Equivalent flat plate area = 0.67 ft 2 

Vertical projected area = 1.35 ft 2 

CG height above waterline = -0.23 ft 
CG fuselage station = 1 .06 ft 

CG position right of buttline = 0.00 ft 

Ixx= 0.719 slug ft 2 
Iyy= 1.433 slug ft 2 
Izz = 2.594 slug ft 2 
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** INPUT DATA CONTINUED (screen 2 of 8) *** 



Main Rotor 

Number of blades = 2 

Rotor radius = 2.83 ft 

Blade twist = 0.00 degs 

Blade airfoil = VR-12 

Blade lift curve slope = 5.73 

Blade weight = 0.42 lbs 

Rotational velocity = 164.00 rads/sec 



Blade grip length = 0.27 ft 

Hinge offset = 0.00 ft 

Flapping moment of inertia = 0.034 slug ft 2 

Hub height above waterline = 0.75 ft 

Hub fuselage station = 1 .25 ft 

Hub position rt of buttline = 0.00 ft 

Mast incidence = 0.00 deg 



*** INPUT DATA CONTINUED (screen 3 of 8) *** 



Tail rotor (zeros if using NOTAR) 
Number of blades = 2.0 

Blade chord = 0.01ft 

Blade radius = 0.52 ft 

Lift curve slope = 5.73 
Rotational velocity = 754.00 rad/sec 



Hub fuselage station = 4.730 ft 

Hub position rt of buttline = -0.188 ft 

Delta-3 angle = 0.00 deg 

Blade twist = 0.00 deg 

Hub height above waterline = 0. 1 56 ft 

Flapping moment of inertia = 0. slug ft 2 



*** INPUT DATA CONTINUED (screen 4 of 8) *** 



NOTAR (zeros if using tail rotor) 
Height above waterline = 0.0 ft 2 

Fuselage station = 0.0 ft 2 

Position right of buttline = 0.0 ft 2 



NOTAR boom diameter = 0.0 ft 2 

Swirl angle at boom = 0.00 deg 
Maximum thruster force = 0.0 lbs 

Thrust fuselage station = 0.0 ft A 2 



*** INPUT DATA CONTINUED (screen 5 of 8) *** 



Wing 

Area = 0.0 ft 2 

Span = 0.0 ft 

CL = 0.00 

CDo = 0.0000 

Tip cord = 0.0 ft 

Root cord = 0.0 ft 

Rotor downwash ratio = 0.00 

Fuselage downwash ratio = 0.00 



Wing efficiency factor = 0.00 
Zero lift angle = 0.00 deg 
Angle of incidence = 0.00 deg 
Lift curve slope = 0.00 
Height above waterline = 0.0 ft 

Fuselage station = 0.0 ft 

Position right of buttline = 0.0 ft 
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*** INPUT DATA CONTINUED (screen 6 of 8) *** 



Horizontal tail 

Area = 0.12 ft 2 

Span = 0.66 ft 

CL = 0.00 

CDo = 0.0050 

Zero lift angle = 0.00 deg 

Angle of incidence = 0.00 deg 



Lift curve slope = 3.4 

Height above waterline = 0.20 ft 

Fuselage station = 3.70 ft 

Position right of buttline = 0.00 ft 

Dynamic pressure ratio = 0.60 
Rotor downwash ratio = 1.50 
Fuselage downwash ratio = 0.06 



*** INPUT DATA CONTINUED (screen 7 of 8) *** 



Vertical tail 
Area = 0.22 ft A 2 

Span = 0.92 ft 

CL = 0.00 
CDo = 0.0050 

Height above waterline = 0.04 ft 



Fuselage station = 4.52 ft 

Position right of buttline = 0.05 ft 

Zero lift angle = 0.00 deg 
Maximum Cl = 0.8 

Dynamic pressure ratio = 0.60 
Lift curve slope = 3.44 



*** INPUT DATA CONTINUED (screen 8 of 8) *** 



Rigging 

Long cyclic pitch/inch defl = 0.03 deg/in 

Lat cyclic pitch/inch defl = 0.01 deg/in 

Collective pitch/inch defl = -0.02 deg/in 

Tail rotor pitch change/defl = -0.09 deg/unit 

Max deflection of control from neutral for NOTAR = 0.00 units 

Displacement of anti-torque control until full rudder = 0.00 units 
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APPENDIX C. MASS-MOMENTS OF INERTIA CALCULATIONS 



This is a listing of the variables used to determine the mass-moments of inertia for the 
Bergen UAV using the technique outlined in Reference 7. 

For Trials #1,2 and 3, the following relation was used: 



/ 



CO 



An 1 



M+S 2 

P M+S 



W z 

[ \Mf_M 

g 



g 



Trial 1 Trial 2 Trial 3 



Component 


Variable 


I xx, Iyy 


Izz 


I xx, Iyy 


Izz 


Ixx, Iyy 


Model Weight 


W M 


15.57 


15.57 


16.44 


16.44 


16.65 


Support Weight 


W s 


2.44 


1.61 


0 


0 


0 


Distance to support CG 


Z s 


6.584 


6.584 


0 


0 


0 


Distance to model CG 


Zm 


13.23 


13.41 


14.01 


13.45 


4.0 


Distance to system CG 


Zm+s 


12.33 


12.77 


14.01 


13.45 


4.0 



* Test set-up on Trial 3 lacked sufficient clearance to check the Izz mode, 
where. 



W 5 z s +w u z u 



The resulting system periods are listed below for each mode; 



Trial 


Ixx 


Iyy 


Izz 


1 


4.0 sec 


4.01 sec 


4.0 sec 


2 


4.15 sec 


4.13 sec 


4.10 sec 


3 


2.24 sec 


2.267 sec 


— 
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APPENDIX D. BEAM DEFLECTION EXPERIMENTAL DATA 



This appendix provides the experimental data and results obtained for the stiffness, El, of 
the main rotor blade using the relation below. 



load, lbs. 


deflection, v 


El (10 6 ) 


1 


1.063 


0.008906 


2 


2.094 


0.009042 


3 


3.03 


0.009373 


. 3.5 


3.44 


0.009632 


0.4 


0.5 


0.007573 


0.9 


1 


0.00852 


1.8 


2 


0.00852 


2.6 


3 


0.008204 


3 


3.5 


0.008114 



EI Avg . = 0.008654 (*10 6 ) 



Beam Deflection vs. Load 
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APPENDIX E. STATE-SPACE MATRICES 



The following is a listing of the output state-space matrices calculated using the JANRAD 
Stability program. 



-0.0091 0 


0.4862 


-32.2000 


-0.0020 




-0.1538 


0 


0 


0 


-1.8778 


0 


0.0002 


0 




0 


1.5383 


0 


0.0037-0.1477 


-0.1972 


0 


0.0005 


0.0624 


0 


.0012 


0 


0 


0.9989 


0 


0 




0 


0 


.0478 


0.0020 0 


-0.1538 


0.0000 


-0.0180 


-0.4896 


32.1632 


.0327 


0.0016 0 


-0.1243 


0 


-0.0101 


-0.3941 


0 


.0104 


0 


0 


0.0000 


0 


0 




1.0000 


0 


0 


0 


-0.1270 


0 


0 


0.0075 


0.0029 


0 


-5.2071 




0.0116 


0.0000 


0 


0 












0 


0.3864 


0 


0 












-0.0047 


0.0304 


0 - 


•0.0040 










B = 


0 


0 


0 


0 












0 


0 


0.0052 • 


-0.0037 












0 


0 


0.0042 


-0.0012 












0 


0 


0 


0 












0 -0.0125 


0 


0.0031 
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APPENDIX F. CENTER OF GRAVITY AND MOMENT OF INERTIA 



The following two EXCEL spreadsheets calculate the changes in CG position, as well as, 
Moments of Inertia with the addition of measurement devices. All inputs are made to 
Sheet 1. 



CenteL °f Gravity Spreadsheet 



This spreadsheet will calculate displacement of the Center of Input a "1" into the 3rd column if 
Gravity with the addition of flight-test instrumentation. the component is installed. 

The new gross weight and center of gravity location **lnput values of X, Y, & Z are in inches 

can be changed in the JANRAD input file (pg.3 of 3). from front-centerline of avionics rack. 



Component 


Weights 


0/1 


Weight 


X 


Y 


Z 


UAV 


16.44 


1 


16.44 


14.25 


0 


12.19 


Gear Mod 


5.64 


0 


0 


12.2 


0 


2.1 


Original Gear 


-0.62 


0 


0 


11 


0 


6.5 


Fuel 


0.855 


1 


0.855 


3.5 


0 


7 


IMU 


2 


0 


0 


18 


0 


2.59 


SONALT TX 


0.19 


0 


0 


10.5 


-5 


-0.5 


SONALTCircuit 


0.42 


0 


0 


19 


-7 


3.5 


Temp SONALT 


1.2 


1 


1.2 


3.5 


0 


5.5 


Freewaves 


0.75 


0 


0 


12 


4.5 


0.33 


Power Panel 


1.66 


0 


0 


18 


6.5 


3.5 


Battery 


2.32 


0 


0 


18 


-4 


4 


Camera 


3 


0 


0 


3 


0 


0 


Autopilot 


0.485 


0 


0 


0 


0 


0 


DGPS Antenna 




0 










DGPS Receiver 




0 










Misc. 




0 










Basic Helo Moments [slug-ft2] 




0.0709 


0.3967 


1.88 



New C.G. Position* 



X- 


1.194 


y- 


0.000 


z- 


2.424 



*WRT reference 
point. 



Gross Weight 18.495 lbs. 
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Moment of Inertia Spreadsheet 



This spreadsheet will automatically calculate the 
new Moments of Inertia of the modified UAV based 
upon inputs on the preceding page. 

Note: Do NOT change inputs to this sheet. 



Basic Helo 


Ixx 

0.0709 


lyy 

0.3967 


Izz 

1.88 


[slug-ft 2 ] 


X- 


1.194 


Modified Helo 


0.086 


0.457 


1.926 


[slug-ft 2 ] 


N *< 

If If 


0.000 

2.424 



Component 


Weight 


mass 


X 


m(y'+z" ) 


Y 


m(x 2 +z 2 ) 


z 


m(x^+y 2 ) 


UAV 


16.44 


0.511 


-1.194 


0.232 


0.000 


0.960 


-0.674 


0.728 


Gear Mod 


0 


0.000 


0.856 


0.000 


0.000 


0.000 


9.416 


0.000 


Original Gear 


0 


0.000 


2.056 


0.000 


0.000 


0.000 


5.016 


0.000 


Fuel 


0.855 


0.027 


9.556 


0.542 


0.000 


2.966 


4.516 


2.424 


IMU 


0 


0.000 


-4.944 


0.000 


0.000 


0.000 


8.926 


0.000 


SON ALT TX 


0 


0.000 


2.556 


0.000 


-5.000 


0.000 


12.016 


0.000 


SONALT 

Circuit 


0 


0.000 


-5.944 


0.000 


-7.000 


0.000 


8.016 


0.000 


Temp SONALT 


1.2 


0.037 


9.556 


1.349 


0.000 


4.752 


6.016 


3.403 


Freewaves 


0 


0.000 


1.056 


0.000 


4.500 


0.000 


11.186 


0.000 


Power Panel 


0 


0.000 


-4.944 


0.000 


6.500 


0.000 


8.016 


0.000 


Battery 


0 


0.000 


-4.944 


0.000 


-4.000 


0.000 


7.516 


0.000 


Camera 


0 


0.000 


10.056 


0.000 


0.000 


0.000 


11.516 


0.000 


Autopilot 


0 


0.000 


13.056 


0.000 


0.000 


0.000 


11.516 


0.000 


DGPS Antenna 


0 


0.000 


13.056 


0.000 


0.000 


0.000 


11.516 


0.000 


DGPS Receiver 


0 


0.000 


13.056 


0.000 


0.000 


0.000 


11.516 


0.000 


Misc. 


0 


0.000 


13.056 


0.000 


0.000 


0.000 


11.516 


0.000 



2.122 8.678 6.556 



Gross Weight 18.495 
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